TESS Discovery of a Super-Earth and Three Sub-Neptunes Hosted by the Bright, Sun-like Star HD 108236

TESS is a spaceborne telescope with four cameras, each with four charge-coupled devices (CCDs) with the primary mission of discovering small planets hosted by bright stars, enabled by its high-precision photometric capability in space (Ricker et al. 2014). The Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) regularly calibrates and reduces TESS data, delivering simple aperture photometry (SAP; Twicken et al. 2010; Morris et al. 2017) light curves and presearch data conditioning (PDC; Smith et al. 2012; Stumpe et al. 2012, 2014) light curves that are corrected for systematics. Then, it searches for periodic transits in the resulting light curves using the Transiting Planet Search (TPS; Jenkins 2002; Jenkins et al. 2017) to search for planets. Unlike the box least-squares (BLS; Kovács et al. 2002), which also searches for transit-like pulse trains while not taking into account the correlation structure of noise, TPS employs a noise-compensating matched filter that jointly characterizes the correlation structure of the observation noise while searching for periodic transits. Finally, it delivers the statistically significant candidates as threshold crossing events (TCEs). As members of the TOI working group, we regularly classify these TCEs as planet candidates and false positives. When vetting TCEs as planet candidates, we use the SPOC validation tests (Twicken et al. 2018; Li et al. 2019), such as:

• 1.  
  the eclipsing binary discrimination test to detect the presence of secondary eclipses and compare the depths of odd and even transits to rule out inconsistencies;
• 2.  
  the centroid offset test to determine whether the centroid of the difference (i.e., out-of-transit minus in-transit) image is statistically consistent with the location of the target star;
• 3.  
  a statistical bootstrap test to estimate the false-positive probability of the TCE when compared to other transit-like features in the light curve; and
• 4.  
  an optical ghost diagnostic test to rule out false-positive hypotheses such as instrumental noise, scattered light, or blended light, based on the correlations between the model transit and light curves derived from the core photometric aperture and a surrounding halo.

HD 108236 was among the list of targets observed by TESS with a cadence of 2 minutes and also included in the TESS Guest Investigator (GI) Cycle I proposal (G011250; PI: Walter, Frederick). It was observed by TESS Camera 2, CCD 2 during Sector 10 (UT 2019 March 26–2019 April 22) and TESS Camera 1, CCD1 during Sector 11 (UT 2019 April 22–2019 May 21). The TESS data were processed by the SPOC pipeline. Then, Sector 10 and 11 TESS data and derived products such as the SAP and PDC light curves, including those of HD 108236, were made public on 2019 June 1 (data release 14) and 2019 June 17 (data release 16), respectively.

The first detection of a TCE consistent with a planetary origin from HD 108236 was obtained in Sector 10 TESS data. The TCE had a period of 14.178 days. However, the light curve also had other transit-like features unrelated to the detected TCE, which promoted HD 108236 to a potentially high-priority, multiplanetary system candidate. Sector 11 TESS data triggered three TCEs, one of which had the same period as that from Sector 10. However, the transits of the other TCEs had inconsistent depths. These initial TCEs from individual sectors were vetted as planet candidates with the expectation that a joint TPS analysis of two sectors of TESS data would resolve the ambiguities on the multiplicity and periods of the planet candidates. The multisector data analysis at the end of Sector 13 resulted in the detection of four TCEs with periods of 14.18, 19.59, 6.20, and 3.80 days and signal-to-noise ratios (S/Ns) of 15.3, 16.2, 11.4, and 8.7, respectively. The PDC light curve of HD 108236 from these two sectors is shown in Figure 2. Subsequently, we released alerts on these four TCEs (i.e., TOI 1233.01, TOI 1233.02, TOI 1233.03, and TOI 1233.04) with planet candidate dispositions on 2019 August 26. For the moment, we will refer to these TCEs that have been vetted as planet candidates using the TOI designations.

Zoom In Zoom Out Reset image size

Time-series photometry of a source is inferred from photoelectrons counted in a grid of pixels on the focal plane. The finite point-spread function (PSF) causes nearby sources to be blended. The focus-limited PSF (FWHM of ∼1 to −2 pixels) and the large pixel size (∼21'') of TESS imply that the resulting time-series photometry of a given target will often have contamination from nearby sources.

Blended light from nearby sources can decrease the depth, δ, of a transit by

$$\begin{eqnarray}&&{\delta }^{{\prime} }=\left(1-\frac{{F}_{{\rm{B}}}}{{F}_{{\rm{T}}}+{F}_{{\rm{B}}}}\right)\delta =(1-D)\delta =\left(1-\frac{f}{1+f}\right)\delta ,\end{eqnarray} \tag{ 1 }$$

where δ' is the diluted transit depth and F_B and F_T are the fluxes of the blended and target source, respectively. Here D is dilution, and f ≡ F_B/F_T is the flux ratio of the blended and target objects. The SPOC pipeline corrects the PDC light curves for this dilution of the transits.

The TESS image of HD 108236 from Sector 10 is shown in Figure 3, along with several archival images of the target, including the Science and Engineering Research Council (SERC) J image taken in 1979, the SERC-I image taken in 1983, and the Anglo-Australian Observatory Second Epoch Survey (AAO-SES) image taken in 1994. The apertures that are used to extract the TESS light curves are also shown for Sectors 10 (red) and 11 (purple). Some of the relatively bright neighbors of HD 108236 are TIC 260647148, TIC 260647113, TIC 260647110, and TIC 260647155, which are 77'', 95'', 108'', and 122'' away and have TESS magnitudes of 13.89, 13.73, 12.94, and 11.67, respectively. Due to the large aperture used to collect light from the bright target HD 108236, the total flux from blended sources is roughly f = 1.2% of the photons coming from HD 108236.

Zoom In Zoom Out Reset image size

Detection of periodic transits in photometric time-series data can be due to any of the following:

• 1.  
  an instrumental (systematic) effect;
• 2.  
  the primary (i.e., brightest) star being eclipsed by a companion star (i.e., eclipsing binary);
• 3.  
  a foreground or background star (i.e., gravitationally not associated with the target) aligned with the target being eclipsed by a stellar companion or transited by a planet;
• 4.  
  the primary or one of the fainter (secondary) stars in a hierarchical multiple star system eclipsing each other or being transited by a planet;
• 5.  
  a nearby star (i.e., gravitationally not associated with the target) being eclipsed by a stellar companion or transited by a planet;
• 6.  
  a star being transited by a planet.

Therefore, we individually considered and ruled out the alternative hypotheses in order to ensure that the planetary classification for the origin of the detected transits was not a false positive.

The first false-positive hypothesis was that the transits could be due to an instrumental effect. The orbital periods of TOI 1233.03 and TOI 1233.04 were close to the multiples of the momentum dump period, which occurred every 3.125 days for Sectors 10 and 11, according to the TESS Data Release Notes.⁴⁶ However, the detected transits did not fall near the momentum dumps. In addition, the transit shapes were inconsistent with that of the typical momentum dump artifact (i.e., sudden drop followed by a gradual rise). The difference images also did not show any evidence of scattered light in the vicinity of HD 108236 during the observations of interest. Furthermore, there were many individual transits detected, which made it extremely unlikely that they were produced by unrelated systematic events. This ruled out the instrumental origin of the detected transits.

The transit model fit performed by the SPOC pipeline on the TESS data indicated that the transit was not grazing and that the depth and shape of the transits were consistent with being of planetary nature. This was also confirmed later with our transit model as discussed in Section 3.9. The SPOC data validation also showed that the apparent positions of the TCEs were all within 1 pixel of HD 108236. Nevertheless, the periodic dimming could be due to any of the sufficiently bright sources in the aperture, since transits or eclipses from nearby or physically associated companion stars could be blending into the aperture. In general, dynamical measurements such as transit timing variations (TTVs) could break this degeneracy. However, the small number of transits and the limited baseline (∼60 days) of the detection data did not yet allow TTVs to be used for vetting.

As a result, follow-up observations were needed to rule out the remaining false-positive hypotheses that the transits are on a target other than the brightest target (i.e., primary). In the remainder of this section, we summarize the data we collected to rule out these false-positive hypotheses.

Upon TESS detection, we obtained reconnaissance spectroscopy follow-up data on HD 108236 using the resources of the SG2 subgroup of TFOP at the Cerro Tololo Inter-American Observatory (CTIO) in Chile, including the Network of Robotic Echelle Spectrographs (NRES) of the Las Cumbres Observatory and the CTIO high-resolution spectrometer (CHIRON), as detailed in Table 4.

3.4.1. LCO/NRES

3.4.2. SMARTS/CHIRON

3.4.3. Ruling Out Aligned Eclipses and Transits

The cross-correlation function and the least-squares deconvolution (LSD) line profile inferred from the reconnaissance spectra rule out a well-separated or even partially blended secondary set of lines, constraining any spatially blended companion with different systemic velocities to be fainter than 5% of the primary at 3σ in the TESS band. This flux ratio is linked to the difference of the magnitudes of the blended source, m_B, and the target source, m_T, as

$$\begin{eqnarray}&&{m}_{{\rm{B}}}-{m}_{{\rm{T}}}=-2.5{\mathrm{log}}_{10}f,\end{eqnarray} \tag{ 2 }$$

which implies that the SG2 data rule out spatially blended sources that have different systemic velocities and that are brighter than TESS magnitude 11.9.

Furthermore, through transit geometry, the undiluted depth, δ ≡ (R_p/R_⋆)², of a full (i.e., nongrazing) transit is linked to full and total transit durations. The total transit duration T_tot is the time interval during which at least some part of the transiting object is occluding the background star, whereas the full transit duration T_full is the time interval during which the transiting object is fully within the stellar disk. Therefore, modeling of the full and total transit durations based on the observed transits allows the estimation of dilution of a transit caused by its neighbors. We inferred the dilution consistent with the observed TESS transits using a methodology similar to that discussed in Section 3.9. The marginal posterior of the dilution requires any blended source to be brighter than TESS magnitude 12.1 at 2σ to produce the observed TESS light curve. Therefore, combined with the constraint from the SG2 data, this rules out the hypothesis that the transits could be produced by a faint foreground or background binary. Furthermore, the fact that there are multiple TCEs on the same target implies that the alignment of unassociated background or foreground eclipses or transits is very unlikely (Lissauer et al. 2012).

The reconnaissance spectroscopy data justified further follow-up of the target to obtain precise radial velocities using the SG4 resources of TFOP (see Table 4).

3.5.1. Magellan II/PFS

3.5.2. Ruling Out Stellar Companions

Table 5 summarizes the radial velocity measurements collected by the SG2 and SG4 subgroups of TFOP. The radial velocities obtained using NRES data are consistent with those from Gaia DR2 (Bailer-Jones et al. 2018; Gaia Collaboration et al. 2018), whereas radial velocities inferred from CHIRON observations have a systematic offset.

Figure 4 shows the radial velocity data from NRES, CHIRON, and PFS after subtracting the mean within each data set. Among the three data sets, the PFS data have the smallest uncertainties (∼1 m s⁻¹). However, they also display variations larger than the uncertainties. This is likely caused by the Doppler shifts due to planets validated in this work.

Zoom In Zoom Out Reset image size

The rms of the radial velocity data from NRES, CHIRON, and PFS is 55, 50, and 3 m s⁻¹, respectively. Using the rms of the PFS radial velocity data, we can place a 3σ upper limit of 1450 M_⊕ on the mass of a companion on a circular orbit around HD 108236 with an orbital period less than 1000 days and an orbital inclination of 90°. Furthermore, assuming circular orbits, the PFS data allow us to rule out stellar masses for the objects that have been observed by TESS to transit HD 108236. This is because the observed rms of the PFS data is much smaller than the expected radial velocity semi-amplitude (∼1 km s⁻¹) from a stellar object having a mass larger than ∼13.6 times the Jovian mass.

We note that we did not use the 12 precise radial velocity measurements from PFS to measure the masses of any of the four planets validated in this work. We leave this to a future work (J. T. Teske et al. 2021, in preparation), where a larger set of precise radial velocity measurements from PFS will be used to accurately measure the masses of the validated planets.

The currently available radial velocity data cannot rule out stellar companions at arbitrary orbital periods, eccentricities, and inclinations. Therefore, a remaining false-positive hypothesis would be a hierarchical system containing planets transiting the primary or the secondary. However, the transiting planets would also have to be giants in this case, in order to compensate for the dilution from the companion star. If more than one such giant planet orbited the companion star, the system would be dynamically unstable. The multiplicity of the transiting objects in the system makes this false-positive hypothesis unlikely. Furthermore, as has been shown in Latham et al. (2011), Lissauer et al. (2012), and Guerrero et al. (2020), it is much less likely for a planet candidate to be a false positive in a multiplanetary system than in a system with a single planet. We therefore discarded this false-positive hypothesis based on the observation of four independent TCEs.

In order to rule out aligned foreground or background stars at close separations, high-resolution images are needed. To obtain high-resolution images in the presence of atmospheric scintillation, we used the speckle imaging technique by taking short exposures of the bright target to factor out the effect of atmospheric turbulence. For this purpose, we used the resources of the SG3 subgroup of TFOP and obtained high-resolution speckle images of HD 108236 with SOAR/HRCam and Gemini/Zorro, as detailed in Table 6.

3.6.1. SOAR/HRCAM

Diffraction-limited resolution was obtained via speckle interferometry by using the High-Resolution Camera (HRCam; Tokovinin et al. 2010; Ziegler et al. 2020) at the 4.1 m SOAR telescope by processing short-exposure images taken with high magnification on UT 2020 January 7. The autocorrelation function and the resulting sensitivity curve are presented in the left panel of Figure 5. A contrast of 5 mag is achieved at a separation of 02.

Zoom In Zoom Out Reset image size

3.6.2. Gemini/Zorro

After ruling out binaries and chance alignments for the target, we then proceeded with ruling out the possibility that the transits detected by TESS could be on nearby stars. HD 108236 is the brightest source within a few arcminutes in its vicinity. Given the depth of the transits observed by TESS (0.302 ± 0.031 ppt,⁴⁹ ${0.517}_{-0.040}^{+0.036}$ ppt, 0.889 ± 0.053 ppt, and 1.175 ± 0.069 ppt), the transit depth would have to be deeper by a certain amount as given by Equations (1) and (2) if the transit was not on HD 108236, but rather on a fainter nearby target. In order to rule out the hypothesis that any of the transits could be on a nearby target, we collected seeing-limited (i.e., with a PSF FWHM of ∼1 arcsecond) photometric time-series data during a predicted transit for each planet candidate (i.e., TOIs 1233.01, 1233.02, 1233.03, and 1233.04) using the resources of the SG1 subgroup of TFOP, including the LCOGT and MEarth telescopes. Table 7 lists these observations. As will be discussed in Section 3.7.4, one of these observations (UT 2020 March 17) resulted in a tentative detection of a transit on target.

3.7.1. LCOGT

3.7.2. MEarth-South

MEarth-South employs an array of eight f/9 40 cm Ritchey-Chrétien telescopes on German equatorial mounts (Irwin et al. 2015). During the data acquisition for this work, only seven of the telescopes were operational. Data were obtained on three nights: UT 2020 March 3 (egress of TOI 1233.01), UT 2020 March 11 (full transits of TOI 1233.02 and TOI 1233.03), and UT 2020 March 17 (full transit of TOI 1233.01). Figure 6 shows the in-focus and defocused fields of the MEarth-South observation on UT 2020 March 17.

Zoom In Zoom Out Reset image size

All observations were conducted using the same observational strategy. Exposure times were 35 s, with six telescopes defocused to half flux diameter of 12 pixels to provide photometry of the target star, and one telescope observing in focus with the target star saturated to provide photometry of any nearby or faint contaminating sources not resolved by the defocused time series. Observations were gathered continuously starting when the target rose above 3 air masses (first observation) or evening twilight (other observations) until morning twilight. Telescope 7 used in the defocused set had a stuck shutter, resulting in smearing of the images during readout, but this did not appear to affect the light curves. The defocused observations were performed with a pixel scale of 084. A photometric aperture with a radius of 17 pixels was used to extract the photometric time series. Data were reduced following standard procedures for MEarth photometry (Irwin et al. 2007).

3.7.3. Ruling Out Nearby Eclipses and Transits

3.7.4. Ground-based Detection of a Transit

A transit of planet d was tentatively detected on UT 2020 March 17 at a 1 m LCOGT-CTIO telescope. The photometric time-series data had a relatively short pre-transit baseline. Therefore, we excluded these observations from the global orbital model in Section 3.9, in order to avoid biasing the fit. However, we fitted the LCOGT-CTIO data separately and inferred a transit duration of 3.8 ± 0.2 hr and a transit depth of 1.1 ± 0.2 ppt, which are consistent with those inferred from the TESS data. The inferred mid-transit time was 2,458,571.3365 ± 0.0035 BJD, indicating a transit arrival 14 minutes late compared to the linear ephemeris model based on the TESS data. The associated light curve is shown in Figure 7.

Zoom In Zoom Out Reset image size

HD 108236 has also been observed by the Wide Angle Search for Planets South (WASP-South) survey (Pollacco et al. 2006) in SAAO, South Africa. WASP-South, an array of eight wide-field cameras, was the southern station of the WASP transit-search project (Pollacco et al. 2006). It observed the field of HD 108236 in 2011 and 2012, when equipped with 200 mm, f/1.8 lenses, and then again in 2013 and 2014, equipped with 85 mm, f/1.2 lenses. It observed on each clear night, over a span of 140 nights in each year, with a typical 10-minute cadence, and accumulated about 58,000 photometric measurements on HD 108236. We searched the data for any rotational modulation using the methods from Maxted et al. (2011). We found no significant periodicity between 1 and 80 days, with a 95% confidence upper limit on the amplitude of 1 mmag. We did not detect any transits in the WASP data, consistent with the expected small transit depths of 0.302 ± 0.031, ${0.517}_{-0.040}^{+0.036}$, 0.889 ± 0.053, and 1.175 ± 0.069 ppt. Planet e had the deepest expected transit; however, its relatively long period likely precluded any detection. The shallow transits of the inner planets also made them undetectable. To determine which region of the parameter space of transiting planets can be ruled out with the WASP data set, we performed injection-recovery tests using allesfitter, which will be introduced in Section 3.9. We injected planets over a grid of periods of 10.1, 15.1, ..., 140.1 days and radii of 8, 8.5, ..., 22 R_⊕. For each planet, we tried to recover the injected signal using Transit Least Squares (TLS; Hippke & Heller 2019). We find that ∼50% of transiting planets with radii 1.3–2 R_J and periods less than 100 days could have been found in the WASP data. The recovery rate drops to ∼20% for planets with radii ∼1 R_J and periods less than 100 days. In contrast, planets much smaller than Jupiter or those on periods longer than 100 days would remain undetected in the WASP data.

Following the vetting of the planet candidates, we modeled the TESS PDC light curve to infer the physical properties of the orbiting planets. In order to model the photometric time-series data, we used allesfitter (Günther & Daylan 2019, 2020). The parameters θ of our forward model M are presented in Table 8. We assumed a transit model with a linear ephemeris. We assumed a generic, eccentric orbit. For limb darkening, we used a transformed basis q₁ and q₂ of the linear u₁ and quadratic u₂ coefficients as (Kipping 2013)

$$\begin{eqnarray}&&{q}_{1}={({u}_{1}+{u}_{2})}^{2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{q}_{2}=0.5\displaystyle \frac{{u}_{1}}{{u}_{1}+{u}_{2}}.\end{eqnarray} \tag{ 4 }$$

Table 8.  Parameters of the Transit Forward Model

Parameter	Explanation	Prior
q1; TESS	First limb-darkening parameter 1	uniform
q2; TESS	Second limb-darkening parameter 2	uniform
logσTESS	Logarithm of the scaling factor for relative flux uncertainties	uniform
$\mathrm{log}{\sigma }_{\mathrm{GP};\mathrm{TESS}}$	Amplitude of the Gaussian process Matérn 3/2 kernel	uniform
$\mathrm{log}{\rho }_{\mathrm{GP};\mathrm{TESS}}$	Timescale of the Gaussian process Matérn 3/2 kernel	uniform
D0;TESS	Dilution of the transit depth due to blended light from neighbors	fixed
Rn/R⋆	Ratio of planet n, Rn, to the radius of the host star, R⋆	uniform
(R⋆ + Rn)/an	Sum of the stellar radius R⋆ and planetary radius Rn	uniform
cosin	Cosine of the orbital inclination, i	uniform
T0;n	Mid-transit time about which the linear ephemeris model pivots, i.e., epoch, in BJD	uniform
Pn	Orbital period of planet n in days	uniform
$\sqrt{{e}_{{\rm{n}}}}\cos {\omega }_{{\rm{d}}}$	Square root of the eccentricity times the cosine of the argument of periastron	uniform
$\sqrt{{e}_{{\rm{n}}}}\sin {\omega }_{{\rm{d}}}$	Square root of the eccentricity times the sine of the argument of periastron	uniform

Download table as:  ASCIITypeset image

We modeled this red noise along with any other stellar variability in the data using a Gaussian process (GP) with a Matérn 3/2 kernel as implemented by celerite (Foreman-Mackey et al. 2017).

When modeling the TESS data, we use the PDC light-curve data product from the SPOC pipeline. We provide the posterior in Table 10 for nuisance parameters, Table 11 for the parameters of planets d and e, and Table 12 for the parameters of planets b and c. Although our nominal results come from allesfitter, we have also repeated the analysis using EXOFASTv2 (Eastman et al. 2019) as a cross-check in order to confirm consistency. EXOFASTv2 has a dynamical prior that avoids orbit crossings and ensures dynamical stability of the analyzed system. A notable result from this analysis were additional constraints on the eccentricities of the planets enabled by the Hill stability prior. The inferred eccentricities were smaller than 0.287, 0.197, 0.164, and 0.149 at a confidence level of 2σ for planets b, c, d, and e, respectively.

Table 10.  Posterior of the Fitting Nuisance Parameters

Parameter	Value	Unit	Fit/Fixed
D0;TESS	0.0		fixed
q1;TESS	${0.23}_{-0.11}^{+0.19}$		fit
q2;TESS	${0.43}_{-0.29}^{+0.36}$		fit
$\mathrm{log}{\sigma }_{\mathrm{TESS}}$	−7.4845 ± 0.0090	log rel. flux.	fit
$\mathrm{log}{\sigma }_{\mathrm{GP};\mathrm{TESS}}$	−8.56 ± 0.13		fit
$\mathrm{log}{\rho }_{\mathrm{GP};\mathrm{TESS}}$	−1.27 ± 0.28		fit

Download table as:  ASCIITypeset image

Table 11.  Posterior of the Fitting Parameters for Planets b and c

Parameter	Value	Unit	Fit/Fixed
Rb/R⋆	0.01638 ± 0.00095		fit
(R⋆ + Rb)/ab	${0.0895}_{-0.0025}^{+0.0028}$		fit
$\cos {i}_{{\rm{b}}}$	${0.037}_{-0.022}^{+0.015}$		fit
T0;b	${2458572.1128}_{-0.0036}^{+0.0031}$	BJD	fit
Pb	${3.79523}_{-0.00044}^{+0.00047}$	d	fit
$\sqrt{{e}_{{\rm{b}}}}\cos {\omega }_{{\rm{b}}}$	−0.00 ± 0.50		fit
$\sqrt{{e}_{{\rm{b}}}}\sin {\omega }_{{\rm{b}}}$	$-{0.03}_{-0.31}^{+0.27}$		fit
Rc/R⋆	${0.02134}_{-0.00083}^{+0.00094}$		fit
(R⋆ + Rc)/ac	${0.0647}_{-0.0019}^{+0.0021}$		fit
$\cos {i}_{{\rm{c}}}$	${0.022}_{-0.014}^{+0.013}$		fit
T0;c	${2458572.3949}_{-0.0020}^{+0.0025}$	BJD	fit
Pc	${6.20370}_{-0.00052}^{+0.00064}$	d	fit
$\sqrt{{e}_{{\rm{c}}}}\cos {\omega }_{{\rm{c}}}$	−0.01 ± 0.49		fit
$\sqrt{{e}_{{\rm{c}}}}\sin {\omega }_{{\rm{c}}}$	$-{0.11}_{-0.29}^{+0.23}$		fit

Download table as:  ASCIITypeset image

Table 12.  Posterior of the Fitting Parameters for Planets d and e

Parameter	Value	Unit	Fit/Fixed
Rd/R⋆	0.02805 ± 0.00095		fit
(R⋆ + Rd)/ad	${0.0375}_{-0.0010}^{+0.0012}$		fit
$\cos {i}_{{\rm{d}}}$	${0.0136}_{-0.0078}^{+0.0065}$		fit
T0;d	${2458571.3368}_{-0.0013}^{+0.0015}$	BJD	fit
Pd	${14.17555}_{-0.0011}^{+0.00099}$	d	fit
$\sqrt{{e}_{{\rm{d}}}}\cos {\omega }_{{\rm{d}}}$	$-{0.03}_{-0.48}^{+0.51}$		fit
$\sqrt{{e}_{{\rm{d}}}}\sin {\omega }_{{\rm{d}}}$	$-{0.04}_{-0.27}^{+0.21}$		fit
Re/R⋆	${0.0323}_{-0.0011}^{+0.0012}$		fit
(R⋆ + Re/ae	${0.03043}_{-0.00089}^{+0.00100}$		fit
$\cos {i}_{{\rm{e}}}$	${0.0118}_{-0.0073}^{+0.0052}$		fit
T0;e	2458586.5677 ± 0.0014	BJD	fit
Pe	${19.5917}_{-0.0020}^{+0.0022}$	d	fit
$\sqrt{{e}_{{\rm{e}}}}\cos {\omega }_{{\rm{e}}}$	${0.01}_{-0.54}^{+0.50}$		fit
$\sqrt{{e}_{{\rm{e}}}}\sin {\omega }_{{\rm{e}}}$	${0.02}_{-0.29}^{+0.23}$		fit

Download table as:  ASCIITypeset image

Table 13.  Posterior of the Inferred Parameters for Planets b and c

Property	Value
R⋆/ab	${0.0881}_{-0.0025}^{+0.0027}$
ab/R⋆	11.35 ± 0.34
Rb/ab	${0.001443}_{-0.000092}^{+0.000100}$
Rb (R⊕)	1.586 ± 0.098
Rb (Rjup)	0.1415 ± 0.0087
ab (R⊙)	10.08 ± 0.36
ab (au)	0.0469 ± 0.0017
ib (deg)	${87.88}_{-0.87}^{+1.3}$
eb	${0.20}_{-0.14}^{+0.30}$
wb (deg)	190 ± 140
btra;b	0.38 ± 0.24
Ttot;b (h)	${2.30}_{-0.11}^{+0.16}$
Tfull;b (h)	${2.20}_{-0.12}^{+0.16}$
ρ⋆;b (cgs)	1.92 ± 0.17
Teq;b (K)	${1099}_{-18}^{+19}$
${\delta }_{\mathrm{tr};{\rm{b}};\mathrm{TESS}}$ (ppt)	0.302 ± 0.031
Pb/Pc	${0.611768}_{-0.000098}^{+0.000092}$
Pb/Pd	0.267731 ± 0.000038
Pb/Pe	0.193716 ± 0.000031
R⋆/ac	${0.0634}_{-0.0018}^{+0.0020}$
ac/R⋆	15.78 ± 0.49
Rc/ ac	${0.001354}_{-0.000067}^{+0.000076}$
Rc (R⊕)	${2.068}_{-0.091}^{+0.10}$
Rc (Rjup)	${0.1845}_{-0.0081}^{+0.0089}$
ac (R⊙)	14.01 ± 0.51
ac (au)	0.0651 ± 0.0024
ic (deg)	${88.72}_{-0.74}^{+0.82}$
ec	${0.18}_{-0.14}^{+0.34}$
wc (deg)	210 ± 120
btra;c	${0.33}_{-0.21}^{+0.25}$
Ttot;c(h)	2.913 ± 0.095
Tfull;c(h)	${2.754}_{-0.094}^{+0.100}$
ρ⋆;c (cgs)	1.93 ± 0.18
Teq;c (K)	${932}_{-16}^{+17}$
${\delta }_{\mathrm{tr};{\rm{c}};\mathrm{TESS}}$ (ppt)	${0.517}_{-0.040}^{+0.036}$
Pc/Pb	${1.63461}_{-0.00025}^{+0.00026}$
Pc/Pd	0.437636 ± 0.000052
Pc/Pe	0.316650 ± 0.000046

Download table as:  ASCIITypeset image

Table 14.  Posterior of the Inferred Parameters for Planets d and e and the Host Star

Property	Value
R⋆/ad	${0.0365}_{-0.0010}^{+0.0011}$
ad/R⋆	${27.39}_{-0.82}^{+0.78}$
Rd/ad	${0.001024}_{-0.000046}^{+0.000048}$
Rd (R⊕)	2.72 ± 0.11
Rd (Rjup)	0.2423 ± 0.0097
ad (R⊙)	24.31 ± 0.87
ad (au)	0.1131 ± 0.0040
id (deg)	${89.22}_{-0.38}^{+0.45}$
ed	${0.17}_{-0.12}^{+0.30}$
wd (deg)	${190}_{-130}^{+140}$
btra;d	${0.35}_{-0.21}^{+0.19}$
Ttot;d (h)	${3.734}_{-0.049}^{+0.066}$
Tfull;d (h)	${3.491}_{-0.057}^{+0.061}$
ρ⋆;d (cgs)	1.93 ± 0.17
Teq;d (K)	${708}_{-12}^{+13}$
${\delta }_{\mathrm{tr};{\rm{d}};\mathrm{TESS}}$ (ppt)	0.889 ± 0.053
Pd/Pb	3.73509 ± 0.00053
Pd/Pc	2.28501 ± 0.00027
Pd/Pe	${0.723548}_{-0.000097}^{+0.000090}$
R⋆/ae	${0.02948}_{-0.00086}^{+0.00097}$
ae/R⋆	${33.9}_{-1.1}^{+1.0}$
Re/ ae	${0.000951}_{-0.000043}^{+0.000049}$
Re (R⊕)	${3.12}_{-0.12}^{+0.13}$
Re (Rjup)	${0.279}_{-0.011}^{+0.012}$
ae (R⊙)	30.1 ± 1.1
ae (au)	0.1400 ± 0.0052
ie (deg)	${89.32}_{-0.30}^{+0.42}$
ee	${0.20}_{-0.13}^{+0.30}$
we (deg)	${170}_{-130}^{+150}$
btra;e	${0.36}_{-0.23}^{+0.20}$
Ttot;e(h)	${4.013}_{-0.057}^{+0.080}$
Tfull;e(h)	${3.712}_{-0.069}^{+0.063}$
ρ⋆;e (cgs)	1.92 ± 0.18
Teq;e (K)	${636}_{-11}^{+12}$
${\delta }_{\mathrm{tr};{\rm{e}};\mathrm{TESS}}$ (ppt)	1.175 ± 0.069
Pe/Pb	5.16220 ± 0.00084
Pe/Pc	3.15806 ± 0.00046
Pe/Pd	${1.38208}_{-0.00017}^{+0.00019}$
Limb-darkening u1;TESS	${0.40}_{-0.24}^{+0.22}$
Limb-darkening u2;TESS	${0.06}_{-0.27}^{+0.36}$
ρ⋆;combined (cgs)	1.93 ± 0.17

Download table as:  ASCIITypeset image

We show in Figure 8 the light curve of each planet folded onto its orbital period and centered at the phase of the primary transit, after masking out the transits of the other planets. Because the orbital period of planet d is close to the orbital period of TESS around Earth (∼13.7 days), a large gap is formed in its phase curve. Figure 9 then shows the individual phase curves, along with the posterior median transit model shown with the blue lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HD 108236 is one of the brightest stars that host four or more planets. As shown in the top row of Figure 12, it is the third-brightest (behind Kepler 444, Campante et al. 2015; HIP 41378, Vanderburg et al. 2016) and the fourth-brightest star (behind Kepler 444, HIP 41378, and Kepler 37; Barclay et al. 2013) in the V and J bands, respectively, that is known to host at least four planets. However, out of these, only Kepler 37 is a Sun-like star, making HD 108236 the brightest Sun-like star in the visual band to harbor at least four transiting planets. This property of HD 108236 makes it an interesting and accessible target from an observational point of view regarding future mass measurements, photometric follow-up, and atmospheric characterization of its transiting planets.

Zoom In Zoom Out Reset image size

The bottom row of Figure 12 also shows the radial velocity semi-amplitude (at fixed planet mass and orbital period) divided by the square root of the host star brightness in the V and J bands, respectively, which are denoted by ${K}_{V}^{{\prime} }$ and ${K}_{J}^{{\prime} }$. The x-axes are normalized so that the top target has the value of 1. Being a Sun-like star, HD 108236 falls to the seventh rank, when the comparison is made in the J band, since low-mass stars generate a larger radial velocity signal for a given companion.

The expected radial velocity semi-amplitudes of the four validated planets based on the predicted masses are in the range of 1.3–2.4 m s⁻¹. Given the brightness of the host star, this implies that the system has good potential for mass measurements in the near future. There are ongoing efforts to measure the masses of all validated transiting planets hosted by HD 108236.

Given the current absence of mass measurements of the planets, we use the probabilistic model of Chen (2017) in order to predict the masses of the validated planets. This model takes into account the measurement, sampling, and intrinsic scatter of known planets in the mass–radius plane. As a result, the large uncertainties of the resulting mass predictions are dominated by this intrinsic system-to-system scatter and not by the posterior radius uncertainties of the planets validated in this work.

The masses of planets b, c, d, and e are predicted as 4 ± 2 M_⊕, 5 ± 3 M_⊕, 8 ± 5 M_⊕, and 10 ± 6 M_⊕, respectively. Hence, planet b is likely a dense, rocky planet, whereas planets c, d, and e are sub-Neptunes with a hydrogen and helium envelope whose radius increases going from planet c to e. Atmospheric escape of volatiles is likely to be strongest for the innermost planet b and should decrease toward the outermost planet e.

Once the radius, mass, and, hence, bulk composition of a planet are determined, the next step in its characterization is the determination of its atmospheric properties. The available data on HD 108236 do not yet allow the atmospheric characterization of its planets. However, sub-Neptunes orbiting HD 108236 are favorable targets for near-future atmospheric characterization as we discuss below.

Given the expected launch of the James Webb Space Telescope (JWST), the transmission spectrum metric (TSM; Kempton et al. 2018),

$$\begin{eqnarray}&&\mathrm{TSM}\propto \displaystyle \frac{{R}_{{\rm{p}}}^{3}{T}_{\mathrm{eq}}}{{M}_{{\rm{p}}}{R}_{\star }^{2}},\end{eqnarray} \tag{ 5 }$$

ranks the relative S/N of different planets assuming observations made with the Near Infrared Imager and Slitless Spectrograph (NIRISS; Maszkiewicz 2017) of JWST, assuming a cloud-free, hydrogen-dominated atmosphere.

The largest uncertainty in predicting the TSMs of the planets orbiting HD 108236 arises from the current unavailability of their mass measurements. We use the predicted masses of planets b, c, d, and e in Equation (5) to obtain preliminary estimates of their TSMs. Based on the brightness of the host star, it is expected that the masses of all validated planets will be measured to better than 40%. Therefore, comparing the TSMs of the validated planets to those of all known sub-Neptunes retrieved from the NASA Exoplanet Archive⁵⁰ Planetary Systems Composite Data with mass measurement uncertainties better than 5σ, we find that the sub-Neptunes HD 108236c, HD 108236d, and HD 108236e fall among the top 20. The super-Earth (planet b) is not included in this TSM ranking, because it is not expected to have a hydrogen-dominated atmosphere. We once again emphasize that these rankings are based on the predicted masses and the actual rankings will depend on the mass measurements of the planets.

The logarithms of the relative TSMs of the planets are plotted against their radii in Figure 13, along with those of the known exoplanets (black points) retrieved from the NASA Exoplanet Archive, where the overall normalizations of the TSMs are arbitrary. We only show those known planets that have a measured mass with an uncertainty better than 40%. The three sub-Neptunes of the HD 108236 system are found to be favorable targets for comparative characterization of sub-Neptune atmospheres.

Zoom In Zoom Out Reset image size

It is worth noting that the TSM ranking of the HD 108236 sub-Neptunes improves with decreasing equilibrium temperature, despite the fact that lowering the temperature acts to reduce the pressure scale height.

As can be seen in Equation (5), the TSM is proportional to the third power of R_p, while inversely proportional to M_p. Although it also scales with M_p, the R_p dependence of M_p is weaker than R_p³. Therefore, the TSM is more sensitive to an increase in planetary radius than a drop in equilibrium temperature. In the HD 108236 system, the radii of planets c, d, and e increase with decreasing equilibrium temperature. As a result, the predicted TSM increases from planet c to e.

Furthermore, although HD 108236 is a relatively bright target, its brightness is below the limiting magnitude of NIRISS/JWST (J magnitude of ∼7; Beichman et al. 2014), making it an appealing transmission spectroscopy target for the instrument.

We also note that planets orbiting HD 108236 span a broad range of radius and equilibrium temperature. Figure 14 shows the distribution of radii and equilibrium temperatures of known planets retrieved from the NASA Exoplanet Archive and those of the planets orbiting HD 108236. The wide range of radii and equilibrium temperatures manifested by the planets allows controlled experiments of how stellar insolation affects the photoevaporation of the volatile envelopes of the orbiting planets by controlling for the stellar type and evolution history (Owen & Campos Estrada 2020).

Zoom In Zoom Out Reset image size

In a multiplanetary system, the displacement from a mean motion resonance (MMR),

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{P^{\prime} }{P}\displaystyle \frac{j-k}{j}-1,\end{eqnarray} \tag{ 6 }$$

of adjacent planet pairs characterizes the proximity of the pair to an MMR, where P' and P are the orbital periods of the outer and inner planets, respectively, j is the nearest integer to the orbital period ratio, and k is the order of the nearest MMR. Proximity to an MMR results in TTVs with a coherence timescale (i.e., superperiod) of P_ttv such that

$$\begin{eqnarray}&&\displaystyle \frac{1}{{P}_{\mathrm{ttv}}}=\left|\displaystyle \frac{j-k}{P}-\displaystyle \frac{j}{P^{\prime} }\right|.\end{eqnarray} \tag{ 7 }$$

The HD 108236 system consists of closely packed planets. However, no pair of the validated planets is on an MMR. The proximities and superperiods of the known adjacent pairs in the HD 108236 system are shown in Table 9.

For the first-order interaction between a pair, where k = 1, the amplitude of the TTVs, V and V', can be estimated using the analytical solution (Lithwick et al. 2012)

$$\begin{eqnarray}&&V=P\displaystyle \frac{\mu ^{\prime} }{\pi {j}^{2/3}{\left(j-1\right)}^{1/3}{\rm{\Delta }}}\left(-f-\displaystyle \frac{3}{2}\displaystyle \frac{{Z}_{\mathrm{free}}^{* }}{{\rm{\Delta }}}\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&V^{\prime} =P^{\prime} \displaystyle \frac{\mu }{\pi j{\rm{\Delta }}}\left(-g+\displaystyle \frac{3}{2}\displaystyle \frac{{Z}_{\mathrm{free}}^{* }}{{\rm{\Delta }}}\right),\end{eqnarray} \tag{ 9 }$$

where f and g are coefficients, μ and μ' are the masses of the planets normalized by that of the host star, and ${Z}_{\mathrm{free}}^{* }$ is the conjugate of the complex sum of eccentricity vectors.

No planet pairs in the HD 108236 system are in or near a strong MMR, precluding the generation of large resonant TTVs. However, nonresonant (chopping) TTVs with small amplitudes induced by synodic interactions are possible. Assuming circular orbits and using the predicted masses yield a TTV of ∼5 minutes for both planets d and e. We also confirmed this analytical prediction using an N-body dynamical simulation (Lissauer et al. 2011) of HD 108236 with a length of 5000 days. We note that the planets could have higher TTVs when the circular orbit assumption is relaxed. Hence, with sufficient transit timing precision, planets d and e are likely to be amenable to mass measurements via TTV observations enabled by long-term transit photometry follow-up (Deck & Agol 2015).

Potential three-body resonances due to a hypothetical planet x.—The orbital gaps between planets b and c and between planets c and d are large enough for low-mass planets to exist on stable orbits, as is common among multiplanetary systems discovered by the Kepler telescope. There are many adjacent pairs in the Kepler data set close to the 3:2 MMR, which invokes the possibility of a missing planet in the apparent 9:4 near-resonant gap between the middle pair of HD 108236. While the parameter space for such missing planets is fairly large, we note that resonant chains of three bodies, as is present in systems like TRAPPIST-1 (Gillon et al. 2017) and Kepler-80 (Xie 2013), could be present in HD 108236 owing to yet-undetected planets. This undetected planet x could have either a period of P_x = 9.364 days, which would satisfy

$$\begin{eqnarray}&&0\approx 2{n}_{{\rm{c}}}-5{n}_{{\rm{x}}}+3{n}_{{\rm{d}}},\end{eqnarray} \tag{ 10 }$$

where n_x is the orbital frequency of the hypothetical planet, or a period of P_x = 9.150 days, which would satisfy

$$\begin{eqnarray}&&0\approx {n}_{{\rm{x}}}-3{n}_{{\rm{d}}}+2{n}_{{\rm{e}}}.\end{eqnarray} \tag{ 11 }$$

The resulting 3:2 resonance between this hypothetical planet x and planet d would result in additional TTVs.

To search for evidence of such an additional planet in the TESS data, we used allesfitter's interface to remove the remnant stellar variability from the PDC light curve using a cubic spline and recursive sigma clipping via wotan (Hippke et al. 2019). Then, we ran a TLS search (Hippke & Heller 2019) on this flattened light curve. We recovered all four transiting planets b, c, d, and e. We also detected several additional periodic transit-like signals above an S/N threshold of 5. The most statistically significant of these detections has an epoch of 2,458,570.6781 BJD, period of 10.9113 days, transit depth of 0.23 ppt, S/N of 8.0, signal detection efficiency (SDE) of 6.9, and false-alarm probability of 0.01. We therefore present this as a tentative fifth planet candidate in the HD 108236 system. However, given the large false-positive probability and its dependence on the detrending method, we concluded that instrumental origin cannot be ruled out for this planet candidate. In particular, the stellar density consistent with the transits of this candidate is 0.4 ± 0.3 g cm⁻³, which is inconsistent with the stellar density (1.9 g cm⁻³) inferred in Section 2. This implies that the candidate is likely due to systematics. Given the larger false-positive probabilities of the other TLS detections (i.e., larger than 0.01), we discarded them as likely due to systematics in the TESS data.

TTV analysis of TESS transits.—In order to infer the TTVs consistent with the TESS data, we performed a light-curve analysis independent of that discussed in Section 3.9 using exoplanet (Foreman-Mackey et al. 2020) by relaxing the assumption of a linear ephemeris. The resulting TTVs are shown in Figure 15. Table 15 also tabulates the mid-transit times of the transits detected in the TESS data. We did not detect any significant TTVs given the temporal baseline and timing precision of the transits observed by TESS. Nevertheless, using these TTVs, we were able to constrain the mass of planet e to be lower than 31 M_⊕ at 2σ via the dynamical simulation, which is consistent with the mass predicted via Chen (2017).

Zoom In Zoom Out Reset image size

Stability.—To further test the dynamical integrity of the system, we conducted N-body integrations using the Mercury Integrator Package (Chambers 1999). Our method is similar to that adopted by Kane (2015, 2019) in the study of compact planetary systems discovered by Kepler. The innermost planet of our system has an orbital period of ∼3.8 days. To ensure perturbative accuracy, we therefore used a conservative time step for the simulations of 0.1 days, which is ∼1/40 of the period of the innermost planet. We ran the simulation for 10⁷ yr, equivalent to ∼10⁹ orbits of the innermost planet. For the masses of planets b, c, d, and e, we assumed fiducial values of 3.5, 4.7, 7.2, and 11.1 M_⊕, respectively. The results of the simulation are represented in Figure 16 by showing the histogram of the eccentricities of the four planets for the entire simulation. The results show that the system is dynamically stable, even considering the nonzero eccentricities for such a compact system. However, there is significant transfer of angular momentum that occurs between the planets with time. The two innermost planets have eccentricities that oscillate between 0 and ∼0.13, which can result in substantial changes in the climate of the atmospheres (Kane & Torres 2017; Way & Georgakarakos 2017), known as Milankovitch cycles (Spiegel et al. 2010). The two outermost planets, d and e, remain near their starting eccentricities and so are largely unperturbed through the orbital evolution.

Zoom In Zoom Out Reset image size